Quadratic Interval Refinement for Real Roots 

John Abbott 
February 25, 2013 

^ : Originally presented as a "Poster" at ISSAC 2006 

*~i ' Abstract 

C^ ■ 

^5 ' We present a new algorithm for refining a real interval containing a single real root: the new method combines 

^^ , characteristics of the classical Bisection algorithm and Newton's Iteration. Our method exhibits quadratic con- 

VO ' vergence when refining isolating intervals of simple roots of polynomials (and other well-behaved functions). We 

assume the use of arbitrary precision rational arithmetic. Unlike Newton's Iteration our method does not need to 
^j , evaluate the derivative. 

• ; 1 Introduction 

^ , 

jrt ' Typically the process of approximating the real roots of a univariate polynomial comprises two phases: root isolation 

rt I where the distinct roots are separated into disjoint intervals, followed by root refinement where the approximations 

I ""^ I . of the roots are improved until they are within specified limits. This paper is concerned with refinement, and 

assumes that isolation has already taken place (see for instance pQ). 
T-H , There are many methods for root refinement. Newton's Iteration is a well-known example, and dates back over 

^ ■ three hundred years. Bisection is another refinement method, slower than Newton's Iteration but more robust. This 

["*"■ ' article presents a variant of Regula Falsi which combines the robustness of Bisection with the rapid convergence of 

t^ , Newton's Iteration. The algorithm is implemented as an integral part of CbQA (from version 4.3 onwards). We 

^^ ' call the new algorithm Quadratic Interval Refinement, or simply QIR. 

. ' Prerequisites for using QIR are knowledge of an initial interval / for which the function has opposite signs at the 

CO ■ two end points (both must be rational), and a procedure which evaluates / at a given rational point with arbitrary 

^^ ' precision. Naturally, the function / must be continuous on /. It is best if the interval / contains a single simple 

CN I root of /; furthermore, if / has a valid Taylor expansion about that single simple root then convergence of QIR 

\ I . will ultimately be quadratic. If / contains several roots then QIR will eventually discard all but one of them as 

^ ' refinement proceeds. 
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1.1 Notation 

For a real number x we write round(x) to mean the integer closest to x; either candidate may be chosen in the case 
of a tie. We use the notation I {a, b) to denote the open interval {x GB. : a < x < b}. 

We abuse terminology harmlessly by taking isolating interval to mean an interval I{a, b) for which /(a) and 
f{b) have opposite sign — we do not require that the interval contain just a single root of /. 

By Newton's Iteration we refer to the root approximation method starting from some given xo and generating 
successive iterates using the formula Xn = x„--i — f{x„-i)/f'{x„-i). This method is sometimes called Newton- 
Raphson Iteration. 

By Bisection we refer to the interval refinement method where the interval I (a, b) is refined to either I {a, {a + 
b)/2) or /((a -I- b)/2, b), the choice depending on the sign of /((o -I- b)/2). This method is also called Binary Chop. 

By Regula Falsi we refer to the interval refinement method where linear interpolation is used to predict the 
position of the root, and then one end of the interval is moved to this predicted value so that the new interval still 
contains a sign change. 

2 Motivation 

Why invent a new method when Newton's Iteration works so well? QIR offers a simple unified approach which can 
refine any isolating interval with a rate of convergence and a computational cost comparable to that of Newton's 
Iteration. Thus QIR has several advantages over Newton's Iteration: 

(a) the criteria guaranteeing convergence are straightforward and easy to verify; 



(b) there is no need to evaluate the derivative of the function whose zeroes we are approximating; 

(c) the successive approximants have simple denominators; 

(d) roots are approximated by intervals containing them, so the accuracy of an approximation is quite explicit. 

A significant difficulty with Newton's Iteration is obtaining a suitable starting value given an isolating interval: 
this is the essence of point (a). It may be necessary to use some other interval refinement method prior to 
using Newton's Iteration. Point (b) poses no problem when approximating roots of explicit polynomials since the 
derivative can readily be evaluated, but in other cases it may not be easy to obtain values of the derivative. Point (c) 
is relevant only when using arbitrary precision rational arithmetic. Point (d) is important when a guarantee of the 
accuracy obtained is desired. With some cunning one can alleviate weaknesses (c) and (d) of Newton's Iteration: 
e.g. replacing approximants by similar values having a simpler denominator of about the right size, and using the 
derivative to estimate the width of an interval containing the root. 

Example QIR can never fail, whereas Newton's Iteration can. The polynomial f{x) = x^ — x-\-Q.7 has a single real 
root, but Newton's Iteration fails to converge when given a starting value xo in the range —0.1 < a;o < 0.1; so 
bad starting points are not so rare. 

Note that one should avoid naive use of Newton's Iteration with exact rational arithmetic because the sizes of 
the numerators and denominators of successive approximants can increase rapidly: typically the fc-th approximant 
will have numerator (or denominator) containing 0(d'°) digits where d is the degree of the polynomial whose root 
is being sought. 

Example The polynomial x"' — 2 contains a real root in the interval 1(1,2). To estimate the root with an error 
less than 2^"'^ starting from the initial approximation xq — 1 requires five Newton iterations, and produces 
an approximant having a denominator with over 500 digits. Alternatively starting from xo = 2 requires seven 
iterations, and the final approximant has a denominator with over 17000 digits. In contrast, QIR obtains an 
interval of width 2~^^ after six iterations, and no value in the computation has more than 50 digits. 



3 Description of the Method 



QIR works by repeatedly narrowing the isolating interval until the width is smaller than the prescribed limit. Each 
individual narrowing step receives an initial interval and a refinement factor A'^ G N by which it must reduce the 
interval width; it uses discretized linear interpolation to guess a narrower interval containing the root. If the guess 
was good, the interval is updated; if not, the interval is unchanged, and the narrowing step returns an indication 
of failure. The refinement factor is increased after every successful narrowing step. 

A narrowing step conceptually divides its initial interval into N consecutive equal width sub-intervals. It then 
uses linear interpolation to locate the sub-interval in which it guesses the root to lie. The guess is tested by 
evaluation at the end points of the chosen sub-interval. If the guess was lucky then the initial interval is replaced 
by the chosen sub-interval; otherwise an indication of failure is returned. The case A'^ = 4 is handled specially: the 
initial interval is always replaced by the correct sub-interval, hut failure is indicated if linear interpolation led to a 
bad guess. 

The refinement factor is varied according to the following rule: 

• after a successful narrowing step, N -h- N"^ 

• after a failed narrowing step, if A > 4 then N <— \/7V 

This strategy is inspired by the knowledge that linear interpolation produces approximations whose errors decrease 
roughly quadratically when sufficiently close to a root of a well-behaved function. 

Obviously, if one of the evaluation points happens to be the exact root ^ then this value is returned as an exact 
root. This is possible only if ^ is rational and with "suitable" denominator. 

3.1 An Illustrative Example 

We present here an example where QIR appears to achieve quadratic convergence immediately (but only due to 
"good fortune"), then later it is forced to reduce the refinement factor back to 4 for a few iterations before true 
quadratic convergence can begin. Note how few evaluations are needed to attain quite high accuracy. 

Example Given the polynomial / = 10 ■ x — 1 the root isolator in CoQA produces the interval 7(0, 2) for the 
positive root. Starting from this interval QIR initially enjoys seven successes, increasing the refinement factor 
to 4^^*, before encountering a number of failures which force the refinement factor back down to 4. Then 
once the interval width has become small enough (less than 10~ , after 24 iterations in Ref inelnterval) 
true quadratic convergence begins. For instance, an interval of width less than lo^^""" is obtained after a 
total of 34 iterations, i.e. at most 68 evaluations of /; a further four iterations give an interval of width less 
than 10-10000. 



4 The Algorithm 

Here we present explicitly the algorithm using a pseudo-language. The actual source code is contained in the CbCbA. 
package RealRoots.cpkg and is publicly available as part of the standard distribution of CbCbA (from version 4.3 
onwards) — see the web site j2j. 

4.1 Main routine: Ref inelnterval 

INPUT: / the function whose zero is being sought, 

I{xio, Xhi) an open interval with rational end points xio < Xhi, and for which f{xio) ■ f{xhi) < 
£ > an upper bound for width of the refined interval. 

OUTPUT: I{(,io,(,hi) an open sub-interval containing a zero of / and having width at most e — the end points 
are both rational. Exceptionally the output may be the exact value of the root ^. 

(1) Initialize refinement factor A'' <— 4, and interval / <— I{xio, xm)- 

(2) While width{I) > e do steps (2.1)-(2.4). 

(2.1) Apply Ref inelntervalByFactor to /, / and A'^ — this will usually modify the interval /. 

(2.2) If an exact root ^ has been found, simply return ^. 

(2.3) If success was reported then replace A' •<— N"^ . 

(2.4) Otherwise failure was reported, so if A^ > 4 replace A' <— y/N . 

(3) Return the interval /. 

4.2 Auxiliary procedure: Ref inelntervalByFactor when A^ > 4 

INPUTS: / the function whose zero we are seeking, 

I — I(xio, Xhi) an open interval with rational end points xio < x^t, and for which f{xio) ■ f{xhi) < 
A' the refinement factor, and we assume A^ > 4 — see >i4.3l for the case A' = 4. 

OUTPUT: success or failure; if success then as a side effect the interval I is replaced by a sub-interval whose 
width is 1/A' times that of I and which contains a zero of /. Exceptionally the output may be exact root ^. 

(1) Conceptually divide the interval I into A'^ consecutive sub-intervals of equal width. Compute the width of each 

sub-interval: w — (xm — xio)/N. Define xo — xio, and Xk = Xk-i + w for fc = 1, 2, . . . , n; so these Xi are the 
end points of the sub-intervals. 

(2) Using linear interpolation predict approximate position of the zero: i. e. determine the closest end point x — Xk 

where k = round(Af • f{xio)/{f{xia) - f{xhi))) 

(3) If fix) — then return exact root x. 

(4) Check whether our prediction was good; i. e. the function / changes sign in one of the sub-intervals having x 
as an end point. 

(5) If f{x) has the same sign as f{xio), check the interval to the right: 

(5.1) If f(x -f li)) = then return exact root x + w. 

(5.2) If f(x + w) has the same sign as f{x) then return failure. 

(5.3) Otherwise replace / •<— I{x, x + w) and return success. 

(6) Otherwise /(£) has the same sign as f(xhi), so check the interval to the left: 

If f(x — w) — then return exact root x — w. 

If f(x — w) has the same sign as f(x) then return failure. 

Otherwise replace / <— I{x — w,x) and return success. 

4.3 Auxiliary procedure: Ref InelntervalByFactor when N = 4 

We give a verbal description rather than pseudo-code for Ref inelntervalByFactor with A'^ = 4. Use linear 
interpolation to predict which of the 5 sub-interval end points is closest to the root. Refine the input interval using 
Bisection twice. If one of the end points of the refined interval agrees with our prediction, return success, otherwise 
return failure. Note that the initial interval is always narrowed by a factor of 4. 



4.4 Termination of the Algorithm 

We show that the algorithm always terminates. In Ref inelnterval if A'^ = 4 then the interval is always narrowed; 
if A'^ > 4 then either the interval is narrowed (by more than a factor of 4), or N is reduced. Since Ref inelnterval 
terminates when the interval has width less than e, and since each narrowing reduces the width by at least a factor 
of 4, there can be only a finite number of narrowings. If no narrowing occurs (i.e. A*' > 4 and failure was reported) 
then A'^ is reduced, but N is never smaller than 4, so there can be only finitely many iterations in which the interval 
is not narrowed. 

4.5 Quadratic Convergence 

The mathematical justification underlying the convergence rate of QIR is very similar to that underlying Newton's 
Iteration. We assume that the isolating interval / contains a single simple root ^ and that the function / admits a 
Taylor expansion centred on ^ valid in an open neighbourhood: i.e. for sufficiently small 5 we have 

f{^ + S) = Ci-S + 02-6" + 0(5''+') 

for some exponent fe > 2 and suitable constants Ci and C2; moreover Ci 7^ since ^ is a simple root. 

Let Xio and Xm be the two end points of the interval / and e be its width. When e is small enough we have: 

fiXlo) = Cl ■ {Xlo - + C2 ■ {Xlo ~ 0" + 0{xio - 0"^^ 
f{Xhi) = Cl■{xu^-S.)+C2■{Xhi~S.f + O{xM-0''^^ 

Estimating ^ by linear interpolation gives ^ = xio + A • [x^i — xio) where A = f{xic)/{f{xio) — f{xhi)). Observe that 
f{xio) ~ fixhi) = Cie + C2?7 + 0(e''+^) where \r]\ = \{xio - i)*' ~ {xm - C)*"! < e*"- Substituting and simplifying 
gives ^ = xio + {£, — xio) — -(^{xio — C)(f ^ i^io — O''^^) + 0{e''^^). Regarding cubic and higher order terms in e 
as negligible, this tells us that |^ — ^| < e • 2C2/C1. 

Now we show that convergence of QIR is eventually quadratic. We say that the isolating interval / of width e is 
narrow enough for / whenever e ■ C2/C1 < 1/8. Observe that, starting from any isolating interval, this condition 
will be satisfied after only finitely many iterations. Henceforth we shall assume / is narrow enough. 

When applying the algorithm Ref inelntervalByFactor to the parameters (/, I, N) we shall say that A'' is small 
enough for / if A'' • £ • C2/C1 < 1/2. Note that A' = 4 is always small enough because we have assumed that / is 
narrow enough. 

Suppose that algorithm Ref inelntervalByFactor is called on {f,I,N) where A'' is small enough for /. The 
condition e ■ C2/C1 < 1/8 implies that linear interpolation produces an estimate £, whose distance from ^ does 
not exceed e^ • 2C2/C1 < e/N . Hence step (2) of Ref inelntervalByFactor makes a good prediction, and so 
success will be returned. Denote the narrowed interval by I' . Now, the next iteration of Ref inelnterval will call 
Ref inelntervalByFactor on the parameters {f, I' , N^). Since I' is narrower than / it is surely narrow enough for 
/. Moreover, one may easily verify that A^^ is small enough for /'. Hence, by induction, all subsequent calls to 
Ref inelntervalByFactor will return success, and convergence is therefore quadratic. 

In contrast, while still assuming that I is narrow enough for /, if Ref inelntervalByFactor is called when A'^ is 
not small enough for / then failure may occur. Since A'^ is not small enough, we must have A'^ > 4. So when failure 
occurs the refinement factor N is reduced. Hence, there can be only finitely many failures before A'' becomes small 
enough for I, and then guaranteed quadratic convergence ensues. 

4.6 Weaknesses 

A weakness of QIR as described here is the cost of evaluating / exactly at various rational points. For instance, if 
/ is a high degree polynomial and the evaluation point a:: is a rational with large denominator then f{x) is likely to 
have an enormous denominator: more or less the n-th power of the denominator of x where n is the degree of /. 

An interesting solution would be to use high precision floating point arithmetic. If variable precision floating 
point numbers are available then it would be sufficient to evaluate / approximately to obtain values with at least 
m + 2 bits of precision where m — log2(A') the logarithm of the refinement factor. If the relative errors in f{xio) 
and f{xhi) are at most 2'™~^ then step (2) in Ref inelntervalByFactor will compute the index k differing by at 
most 1 from the value one would have obtained using exact rational arithmetic. 

As presented here QIR may calculate an interval much narrower than the specified limit. This is a (minor) 
defect because it means the last two evaluations of / will be at rational numbers with needlessly large denominators 
(with up to twice as many digits as truly required), and thus costing more than is strictly necessary. In practice 
one can simply reduce the value of the refinement factor N in the final iteration of Ref inelnterval to avoid the 
overshoot. 



5 Experimental Results 

We present here a small selection of experimental results. The computer used was a Macintosh G5 (running 
MacOS X 10.4) with 2GHz processors, and 3.5G bytes of RAM. The software used was CbCbA 4.6 which relies 
upon the excellent library GMP (version 4.2, see i3^) for arbitrary precision arithmetic. All timings are measured 
in seconds. 

The real root isolation code in CbCbA. is a simplistic implementation of the method described in pp, which is 
based on Descartes's Rule of Signs. Their implementation appears to be faster than ours in OoQA. 

We describe briefly the polynomials used for the timing tests. The first two polynomials are adapted from 
the FRISCO test suite (see ^ and [5]) since CoCoA cannot represent complex numbers. Both polynomials are 
designed to be challenging for a root isolator, and indeed the real root isolator implemented in OaOcA fails with 
"Recursion too deep". However 7(0, 1) is a valid isolating interval for both polynomials, and it was used as input 
to QIR. The first polynomial is /i = {{c^x^ - 3)" + c*x^^){c'^x^ - 3) for c = 10^""; there are four complex roots 
quite close to the real positive root, and this delays the onset of quadratic convergence. The second polynomial is 
/2 = a;^° + (lO^^x — 1)"^ and has two complex roots very close to the real positive root again delaying the onset of 
quadratic convergence: note how quickly 10000 digits are obtained compared to the time needed for 1000 digits. 

The third and fourth polynomials have a similar structure; all their roots are real by construction, /s has degree 
32 and is the minimal polynomial of \/T7E + \/190 + \/T95 + \/398 + \/1482, while /4 has degree 128 and is the 
minimal polynomial of ^99 + \/627 + \/66T + VffS + ^929 + V1366 + V1992. These polynomials are interesting 
because they both have a root which is "almost integer" : /a has a root close to 45, and /i has a root close to 10. 

In the table below, the column headed "#roots" indicates how many intervals were refined, the column headed 
"isolate" is the time taken by the root isolator, and the other columns give the times taken to refine all the intervals 
found by the isolator to a width not greater than e. 
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The GMP library offers functions for computing the integer part of roots of big integers {viz. mpz_sqrt and 
mpz_root); these functions use Newton's Iteration specialized to the computation of fc-th roots. We compared QIR 
with these functions in the following way. To compute d digits of \/5 using mpz_sqrt we passed as argument the 
value 5 x lO^'*. Similarly to compute d digits of v^ and \/2 we called mpz_root with arguments 3 x lO^'' and 
2 X 10^"^ respectively. The timings for QIR are for refining the root interval to a width not exceeding lO""* — the 
initial isolating interval was [0,4] in both cases. The square root function mpz_sqrt is especially efficient: e.g. it 
computed lO'^ digits of \/5 in 8.2 seconds while QIR needed 35 seconds to achieve the same accuracy. In contrast, 
we see that QIR is quite competitive with mpz_root; here are the results: 
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6 Conclusions 

QIR is a simple algorithm combining the speed of Newton's Iteration with the robustness of Bisection. It is easy to 
implement, and the implementation supplied with OaOaA demonstrates its genuine competitiveness with Newton's 
Iteration (except perhaps for computing square roots). 
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